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We examine the steady state properties of binary systems of driven inelastic hard spheres. The 
spheres, which move under the influence of gravity, are contained in a vertical cylinder with a vibrat- 
ing base. We computed the trajectories of the spheres using an event-driven molecular dynamics 
algorithm. In the first part of the study, we chose simulation parameters that match those of ex- 
periments performed by Wildman and Parker yj . Various properties computed from the simulation 
including the density profile, granular temperature and circulation pattern are in good qualitative 
agreement with the experiments. We then studied the effect of varying the mass ratio and the 
size ratio independently while holding the other parameters constant. The mass and size ratio are 
shown to affect the distribution of the energy. The changes in the energy distributions affect the 
packing fraction and temperature of each component. The temperature of the heavier component 
has a non- linear dependence on the mass of the lighter component, while the temperature of the 
lighter component is approximately proportional to its mass. The temperature of both components 
is inversely dependent on the size of the smaller component. 

PACS numbers: 45.70.Mg, 47.20.Bp, 47.27.Te, 81.05.Rm 



I. INTRODUCTION 

Granular systems exhibit many properties that are dif- 
ferent from systems composed of elastic particles. For 
example, driven granular systems display standing and 
travelingwaves 0, 0, , oscillons [f| , heaping, and con- 
vection |y, ML In addition, granular mixtures show size 
segregation Q and steady-state kinetic energies that are 
not equal for each component in the mixture @- This 
departure from equipartition is not unexpected, but it 
is one of the most striking differences between elastic 
and inelastic systems. Understanding the properties of 
mixtures is particularly important for granular systems 
since, unlike molecular systems, they are never com- 
pletely monodisperse. 

Theoretical studies of granular systems have focused 
on two distinct classes. One consists of systems that 
are not driven or heated. The initial energy decays over 
time as a result of inelastic collisions. During this "cool- 
ing" process there is a period during which the density 
is homogeneous. Several workers have presented kinetic 
theories [lJJ, EH jl^l and mean-field theories based on 
Maxwell models [K| to describe the properties of mixed 
granular systems in this homogenous cooling state. In 
the other class of systems, an energy source, such as a vi- 
brating wall, is present. This leads to a non-equilibrium 
steady state that has been studied by several researchers 
[lot \\AY Ha . Ha ] . In both classes the components are pre- 
dicted to have different kinetic energies, or granular tem- 
peratures, that depend on the mass, size, and restitution 
coefficient of the grains 0, 0] . 
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Recently, two-dimensional |17I ] and three-dimensional 
systems Q, [t| of driven, granular mixtures have been 
studied experimentally. Losert et at, who first reported 
the difference in granular temperature |9j , observed that 
the velocity distributions deviated from those observed in 
systems with elastic collisions. Feitosa and Menon stud- 
ied density distributions and granular temperature pro- 
files in two-dimensional systems with and without gravity 
Wildman and Parker have studied the convection 
patterns, density distributions, and temperature profiles 
in three-dimensional systems lj . These studies deter- 
mined that the heavier particles are at a higher granular 
temperature than the lighter particles. In both two and 
three-dimensional systems, the ratio of the temperatures 
varies as the relative proportion of the heavy and light 
particles is changed. The temperature ratio, however, is 
independent of the inelasticity of the particles 0, 0] . 

While the experimental techniques employed in the 
studies cited above have provided many useful insights 
into granular behavior, they cannot easily isolate the ef- 
fects of particle mass, size, and inelasticity. Theoretical 
and computational methods are useful in this respect. 
Molecular dynamics simulations of granular mixtures can 
accurately reproduce the phenomena observed in experi- 
ment |l8l |19| , while providing information on the effects 
of the individual properties mentioned above. For exam- 
ple, Paolotti et al. jnj and Barrat and Trizac [2(j inves- 
tigated the effects of rotation, mass ratio, and relative 
density in two-dimensional vibrated systems. Mixtures 
have also been studied under two-dimensional shear flow 
conditions [21]]. Simulations of the homogeneous cool- 
ing state in two-dimensional systems are consistent with 
experiment and theory [2j,|20, 21]. 

It is important to stress that the conclusions drawn 
from simulations of two-dimensional systems cannot nec- 
essarily be extended to three-dimensions. In particular, 
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the system boundaries have a much larger influence in 
three-dimensions, as recently demonstrated by Talbot 
and Viot [3. 

This paper presents a three-dimensional, event-driven, 
molecular dynamics simulation of a mixture of inelastic 
hard-spheres. The simulation methodology is discussed 
in Section [H] followed by a comparison to the available 
experimental results, Section IIII Al Finally, the effects 

and 



of isolated changes in the mass ratio (Section If 11 B|) 
size ratio fSection 1111 C|) on the energy distribution and 
component temperatures are examined. 



II. MODEL AND SIMULATION 

The three-dimensional system ^| (Figure ^) contains 
a mixture of inelastic hard spheres in an infinitely tall 
cylinder of radius R under the influence of gravity. The 
mixture is composed of ri\ spheres of mass m\ and di- 
ameter d\ and 7J2 spheres of mass m-2 and diameter di. 
Energy is injected into the system by means of the base 
of the cylinder, which vibrates in a symmetric saw-tooth 
waveform of amplitude A and frequency v. 

Three kinds of collisions occur in the system: particle- 
particle, particle-wall and particle-base. The post- 
collisional velocities (v^ ; and j) resulting from a colli- 
sion between a particle from component a and a particle 
from component [3 with masses m a and mp and initial 
velocities v tti i and v^j, respectively, are given by 
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The unit vector between the centers of the colliding parti- 
cles is n, and c is the appropriate restitution coefficient. 
While Equations ^j] and II bl conserve momentum, they 
imply an energy loss of 

A£=-^(l-c 2 )[(v a!i -v w ).nf. (2) 

where /i = m a mp/ (m a + mp) is the reduced mass of the 
particles involved in the collision. 

In a binary mixture it is generally necessary to spec- 
ify three restitution coefficients for particle-particle col- 
lisions: en and C22 for intra-componcnt collisions and 
Ci2(= C21) for inter-component collisions. Several au- 
thors have reported on random [H IH and velocity- 
dependent |23l l24j restitution coefficients. Luding and 
McNamara have proposed a contact time model that 
also leads to a variable restitution coefficient [25]. For 
simplicity, we have chosen to use a constant value for 
the restitution coefficients. We further assume that 

Cll = C 22 = C12 = c. 

Particle-wall collisions are governed by 
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FIG. 1: The three-dimensional system consisting of an in- 
finitely tall cylinder of radius R, and a mixture of hard spheres 
with sizes di and di and masses mi and m2. The base of the 
cylinder is shaken with a symmetric saw-tooth waveform with 
amplitude A and frequency v. 



where c a . w is the appropriate restitution coefficient for 
component a, and f is the radial unit vector. We assume 
that the restitution coefficient for collision with the wall is 
constant for both species and c± tW = C2. w — c w . Particle- 
base collisions are governed by 

v a,i = y a ,i - (1 + C a ,b) [(v Q ,i - V w ) • k] k (4) 

where c a ,b is the appropriate restitution coefficient for 
component a, v w is the velocity of the base at the instant 
of collision, and k is the unit vector in the z-direction. 
The restitution coefficients for collisions with the base are 
also assumed to be constant and equal for both species 
(i.e. c li6 = c 2ib = c fc ). 

A phenomenon similar to inelastic collapse can be ob- 
served in these simulations. For certain ranges of the ve- 



3 



locity a given particle will collide repeatedly with the side 
wall. As its energy is dissipated, the particle approaches 
the wall ever more closely. This is accompanied by an 
increase in collision frequency that eventually "freezes" 
the simulation. To prevent this phenomenon from occur- 
ring, a small impulse is imparted to the particle toward 
the center of the cylinder once its radial velocity falls 
below a certain value. This method has been used previ- 
ously [3, and the threshold value was set such that the 
injected energy does not discernibly influence the simu- 
lation output. It is also possible for a particle to come 
to rest on the base for a time corresponding to 1/2 a cy- 
cle of the vibration. To avoid this possibility, the sign 
of the z-component of the velocity for the particle is in- 
verted when the velocity of the colliding particle is found 
to match the velocity of the base. This condition was 
found to occur once every 12.5 x 10 6 collisions for the 
parameter values used in this paper. Thus, this method 
causes little perturbation in the simulation output. 

We calculated a number of properties from the particle 
positions and velocities generated by the simulation. The 
packing fraction r\ a for component a is defined as 
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where n a is the number of particles of component a in 
the volume element V, and v a = 7rd^/6 is the volume of 
a particle of this component. Another property that we 
calculated is the kinetic energy or granular temperature, 
T a , of each component using the following equations. 



T a = 



T a , y = m a (v a y ) 
T a , z = m a {v 2 a z ) 



(6a) 
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The "partial" temperatures in the x-,y-, and z-directions 
are T a ^ x , T a ^ yi T az , and the angular brackets denote a 
time average over all particles of component a. 

Our first objective was to model the experimental sys- 
tem studied by Wildman and Parker Q . We chose simu- 
lation parameters that correspond to those of the exper- 
iment, i.e., a cylinder of diameter 145mm that is shaken 
at 50 Hz with an amplitude of 1.74 mm. The acceleration 
due to gravity is taken as g = 9.81m/s 2 . The restitution 
coefficient for particle-particle collisions is c = 0.91, for 
particle-wall collisions is c w — 0.68, and for particle-base 
collisions is c b — 0.88. 

We performed simulations in which we varied the rela- 
tive proportions of large and small particles while main- 
taining enough particles to cover the base of the cylin- 
der with a monolayer for comparison to the experimental 
work of Wildman and Parker [lj . We then performed ad- 
ditional simulations to examine the effect of varying the 



mass ratio 7712/ mi and the size ratio c^/^i independently. 
In the discussion of these three studies, component 2 will 
always refer to the smaller and/or lighter component. 

The particles of each component were randomly placed 
in the cylinder with random velocities. We then equili- 
brated all systems for approximately 5000 collisions per 
particle. Data were collected over 2.4 x 10 4 collisions per 
particle at intervals of approximately 10 collisions per 
particle. These data were then averaged for each compo- 
nent to obtain a representation of the system at steady 
state. 



III. RESULTS AND DISCUSSION 
A. Systems with varying composition 

For the simulations discussed in this section, the size 
ratio of the two components was set at dij d\ = 0.8 and 
the mass ratio was set at TO2/TO1 = 0.512. Systems with 
the following compositions were then studied: n\ = 525, 
n 2 = 270; m = 350, n 2 = 540; and m = 175, n 2 = 810. 



1. Velocity Field 

The velocity fields of components 1 and 2 for a system 
with n\ — 525 and n% — 270 particles are shown in Fig- 
ures and [2b, respectively. It can be seen that both 
components circulate in a pattern that rises in the cen- 
ter of the cylinder and falls at the walls. This convection 
pattern has been observed previously in both experiment 
and simulation 0, 0, [2(j . The patterns shown here are 
very similar to those reported by Wildman and Parker p| 
with the center of the convection for both species present 
at a radius of approximately 50 mm and a height of ap- 
proximately 40 mm. 



2. Packing Fraction 

The packing fraction, as a function of radius and height 
for each component, is shown in Figure [3J It can be seen 
that there is a density gradient in both the radial and 
vertical directions for both components, just as is ob- 
served in monodisperse systems |0, |2^ . The data from 
the simulation are qualitatively similar to those observed 
in experiment except that the simulation shows a higher 
concentration at the bottom of the cylinder near the wall 
for both components. It should be noted that while the 
maximum density occurs at the same point for both com- 
ponents, component 2 obtains a greater height than com- 
ponent 1. This trend is opposite that observed for size 
segregation in weakly tapped systems where the larger 
particles rise above the smaller ones ||. 

The radially averaged packing fractions as a function 
of height are shown in Figure 0] The three curves corre- 
spond to the three relative fractions of component 1 and 
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FIG. 2: Velocity field of (a) component 1 and (b) component 2. Simulation conditions as follow: ni 
dijd\ = 0.8, m^/mi = 0.512, all other conditions as given in the text. 
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FIG. 3: Contour plots of the packing fraction of (a) component 1 and (b) component 2. Contours correspond to a packing 
fraction change of 0.002 in (a) and 0.0011 in (b). Simulation conditions as given in Figure|5| 
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FIG. 4: Packing fractions of (a) component 1 and (b) com- 
ponent 2 for the compositions: I tii = 525, ni = 270; © 
ni = 350, ri2 = 540; ▼ ni = 175, n2 = 810. Other simulation 
conditions as given in Figure 

component 2. In all cases, the packing fraction of both 
components increases steeply at small heights, reaches 
a maximum and decays at large heights. The packing 
fraction of component 1 decreases and it increases for 
component 2 as the relative amount of each component 
is changed. The changes in the relative fractions only 
affects the magnitudes of the packing fractions. There is 
no noticeable variation in the details of the packing frac- 
tion profiles (i.e. the position of the maximum, the rate 
of decay, etc.) as the relative amounts of the components 
are changed. 



3. Temperature 

The granular temperature of each component is also 
studied as a function of the relative proportions. Fig- 
ure [S] shows contour plots of the granular temperature of 



component 1 in the x- (Figure , y- (Figure [SJ}) , and 
z- (Figure Et) directions for the system presented in Fig- 
ure [3 The symmetry evident in the x- and y-directions 
is produced by the unbiased introduction of energy into 
these directions by particle-particle collisions. These two 
partial temperatures decay rapidly in both the radial and 
vertical dimensions from a maximum near the center of 
the cylinder, close to the vibrating base. The z-direction, 
however, is different because of the bias introduced by the 
vibrating base. This partial temperature decays in the 
vertical dimension with very little radial dependence at 
small heights. At larger heights, the center of the cylinder 
is slightly warmer than the surrounding area, as would 
be expected given the toroidal flow profile presented in 
Figure |21 

Figure [S] shows the height dependence of the radially 
averaged, partial granular temperatures of each compo- 
nent. It can be seen that the height profiles of the tem- 
perature are similar for both components at all relative 
proportions. The temperatures in the x- and y-directions 
show the formation of a maximum for both components 
as height increases. These maxima occur very close to 
the height that corresponds to the maximum in packing 
fraction. The increase in the partial temperatures can 
then be attributed to an increase of particle-particle col- 
lisions that inject energy into the x- and y-directions, in- 
creasing the corresponding temperatures. Figure also 
shows that the temperature in the z-direction is larger 
than that in the other two directions. The minimum ob- 
served as height increases has been predicted by Brey 
et al. for systems in which the particles do not interact 
with a top barrier, but are under the influence of gravity 
|27|. Wildman et al. observed the minimum the experi- 
mental systems [2(| [28| , and Ramirez and Soto presented 
a hydrodynamic theory that addresses this phenomenon 

Figures & and[|J> also show that the temperatures of 
the two components change as the relative proportion of 
the two components changes. The temperature in the 
three spatial directions decreases as the relative propor- 
tion of the larger particles is decreased. As the fraction 
of component 1 decreases, the number of collisions with 
the smaller particles increases, causing the temperature 
of the larger particles to decrease. Also, the temperature 
of the entire system decreases because component 2 does 
not gain as much kinetic energy from the base since it is 
lighter than component 1. The decrease in the temper- 
ature, however, does not change the height at which the 
extrema in T x , T y and T z are observed. Only the magni- 
tude of the measured temperature appears to change in 
these systems (Figure andEJa). 

While the trends observed in the temperature of the 
two components are similar, there are differences between 
the two components. The temperatures of component 1 
(FigureHOi) are greater than those of component 2 (Figure 
in all the systems, an effect particularly pronounced 
in the z-direction. The temperatures in the x- and y- 
directions differ only slightly. It is difficult to determine 
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FIG. 5: Contour plots of the temperature in the (a) x-direction, (b) y-direction, and (c) z-direction for component 1. Contours 
correspond to a change of 11.5 fj,J. Simulation conditions as given in Figurc|5| 



whether the differences in temperature between the two 
components are dominated by the differences in size or 
mass from these data. Thus, we conducted further sim- 
ulations to determine the individual effects of mass and 
size. 



B. Systems with varying mass ratio 

To determine the effects of particle mass, we simulated 
systems with mass ratios m^jriix = 0.01, 0.125, 0.25, 
0.5, and 1.0 at constant size, dijd\ = 1.0, and relative 
fraction ria/^i = 1.0 and with a total n\ + rt2 = 1050 
particles in the system. 



1. Energy Distribution 

We expect changes in the mass ratio to affect the ex- 
change of energy, and hence the temperature, of each 
component. We therefore computed the average change 
in energy, (AE a )p, of a particle of component a result- 
ing from collisions with particles of component /?. We 
examined the energy exchanges resulting from particle- 
particle, particle-wall and particle-base collisions. The 
data we collected from the simulations for each kind of 
collision are presented in Figured 

As expected, both components lose energy on collision 



with the wall, while they experience a net energy gain on 
collision with the base JTJi). Both components also lose 
energy from intra-component, particle-particle collisions 
(squares in 03). The inter-component, particle-particle 
collisions (circles in 0d) , however, show trends that are 
not intuitively obvious. Component 1 shows a loss of 
energy for all mass ratios, while component 2 shows that 
there may be a loss or gain of energy depending on the 
mass ratio. We obtained a theoretical estimate for this 
quantity by assuming each component has a Maxwell- 
Boltzmann velocity distribution but with a temperature 
specific to the component. The details of the derivation 
are presented in the Appendix. The average energy loss 
for a particle of component a resulting from collisions 
with particles of component (3 and average component 
temperatures (kinetic energies) of T a and Tp is 

(AE a ) p = k B fi (1 + c) (1 + c) '- - 2 

where ks is the Boltzmann constant. The total average 
energy loss per collision between particles of components 
a and [3 with temperatures of T a and Tp is 

(AE a )p + (AEp) a = -k B (1 - c 2 ) T Z a t^ mf3 (8) 

Equation [S] indicates that for inelastic collisions, c < 1, 
the total energy of the colliding pair always decreases, 
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FIG. 6: Temperatures in the M x-, © y-, and A z-directions 
of (a) component 1 and (b) component 2. The three spacial 
directions are shaded the same for each system. Simulation 
conditions as follows: n\ = 525, ri2 = 270 (top curves); n\ = 
350, n2 = 540 (middle curves); ni = 175, 112 = 810 (bottom 
curves) . 



even if the temperatures of the two components are differ- 
ent. For equal temperatures, Equation [7| shows that this 
is also true for the individual energies of the components. 
If, however, the temperatures are different, it is possible 
for the energy of the lighter component to increase on 
average due to collisions with the heavier component. 

We calculated values of (AE a )p for inter-component 
and intra-component, particle-particle collisions for each 
component using Equation The masses of the par- 
ticles and the restitution coefficient are set by the in- 
put parameters, but the temperatures of each compo- 
nent are not known a priori. Therefore, the temperature 
was obtained from the simulation output for each mass 
ratio. The results from Equation {7\ are presented in Fig- 
ure EJa along with the results obtained from the simu- 
lation. It can be seen that the average energy changes 
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FIG. 7: Average energy change per particle per collision for 
(a) particle-boundary collisions and (b) particle-particle col- 
lisions as the mass ratio is changed. Component 1 is repre- 
sented by the closed symbols and component 2 is represented 
by the open symbols. The different kinds of collisions are: M 
intra-component, particle-particle ((AEi)i and (AE^^), ^ 
inter-component, particle-particle ((AEi)? and (A-E/2}i), ^ 
particle-wall, and A particle-base. The solid (dashed) lines 
are the energy loss for component 1 (component 2) as cal- 
culated by Equation Q f° r intra-component (same shade as 
the squares) and inter-component (same shade as the circles), 
particle-particle collisions. 



calculated for the two kinds of particle-particle collisions 
compare favorably with those obtained from the simu- 
lation. There is better agreement for intra-component 
collisions because there is no temperature or mass dif- 
ference between the colliding particles. The predicted 
energy changes for inter-component collisions show the 
same trend that is observed in the simulation results, 
but the magnitude of the change is incorrect. Specifically, 
the equation overestimates the change in energy resulting 
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from inter-component, particle-particle collisions. This 
error probably arises from the inhomogeneities in the 
particle density and temperature (see Figures and EJ , 
which are not accounted for in the Maxwell-Boltzmann 
distribution. 

More generally, we note that as the mass ratio de- 
creases, the energy change associated with any collision 
also decreases. This variation in the energy change will 
affect the bulk properties observed for these systems, 
such as the packing fraction and the temperature dis- 
cussed below. 



2. Packing Fraction 

The packing fraction of each component as a function 
of height is presented in Figure 03 The general behavior 
is the same for both components, and is similar to that 
already discussed in Section lTlI A 21 The packing fraction 
of component 1 varies little with mass ratio changes. In 
particular, the maximum density is at the same height for 
all the systems examined. At greater heights, however, 
component 1 condenses and then expands as the mass 
of component 2 decreases. This phenomenon is more 
clearly visible in the insert of Figure which shows a 
linear relationship between ln(ry) and the height. At large 
altitudes, the slope of the lines in the insert increase as 
component 1 condenses and decrease as component 1 ex- 
pands in the system. This behavior corresponds to the in- 
creased energy loss that is observed for inter-component 
particle-particle collisions (see Figure 0). 

The packing fraction profiles of component 2 undergo 
a much more significant change as the mass ratio is de- 
creased. Specifically, Figure IHJd indicates a steady de- 
pletion of this component from around the maximum, 
with a compensating increase at large altitudes, as the 
particles expand into the upper reaches of the cylinder. 
For mass ratios of 0.25 and below, two local maxima are 
present. These are most distinct for the systems in which 
component 2 has a net gain in energy due to collisions 
with particles of component 1 (mass ratios of 0.125 and 
0.01). Thus, the two maxima are formed as the particles 
of component 2 try to separate from component 1 . The 
increase in the energy forces the particles of component 
2 toward the base and toward higher altitudes. This is 
what is observed in the packing fractions shown in Fig- 
ure[SjD with one maximum very close to the base, and one 
maximum that increases in altitude as the mass ratio de- 
creases. The plots of ln(?7), shown in the insert, display a 
steady decrease in the slope as the mass ratio decreases. 



3. Temperature 

The effect of mass ratio on the temperature in the z- 
direction for components 1 and 2 is shown in Figures 
and respectively. The minimum in the temperature 
is obvious for both particles. It is also easily seen that 
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FIG. 8: Packing fraction of (a) component 1 and (b) compo- 
nent 2 for mass ratios of M ma/mi = 1.0, 9 ira/mi = 0.5, 
♦ milm\ — 0.25, A m^/mi — 0.125, and ▼ milm\ = 0.01. 
All other simulation parameters are as given in the text. The 
inset shows the natural log of the trailing edge of the packing 
fraction for each system. 



the temperature in the z direction of component 1 goes 
through a minimum as the mass ratio of the two com- 
ponents decreases. The temperatures in the x- and y- 
directions (not shown) also follow the same trend. This 
indicates that total temperature for component 1 goes 
through a minimum as the mass ratio is decreased. The 
changes in the temperature coincide with the changes 
observed in the packing fraction (Figure and the en- 
ergy changes for the different kinds of collisions (Fig- 
ure [7| ) . This implies that the changes in the velocities 
of the particles of component 1 affect the temperature, 
just as expected from Equations As shown by Brey 
et al. |27j and Warr et al. [30| . a relationship exists 
between the temperature and the packing fraction in a 
single component system. Extending their results to a 
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FIG. 9: The temperature in the z-direction of (a) component 
1 and (b) component 2 for mass ratios of B mijm\ = 1.0, 
B milm\ = 0.5, ♦ milrax = 0.25, ^ mz/mi = 0.125, and 
▼ mijm\ = 0.01. The temperature calculated from Equation 
[5|is shown as a line corresponding to the data points of the 
same shade. 



multi-component system, we obtain 



rfln(rfo) 

dz 



m a g 
ksT a 



(9) 



This equation holds for high altitudes and restitution co- 
efficients close to 1. Thus, a limiting temperature can be 
calculated for each system using the data presented in 
the insert of Figure^. The results, plotted as horizontal 
lines in Figure 05Jd, correspond well with the asymptotic 
temperatures. This indicates that the decaying edge of 
the packing fraction is a good indicator of the tempera- 
ture at those altitudes. 

We observed trends for component 2 that are very dif- 
ferent from those just discussed for component 1. Figure 
[SJj shows that the temperature decreases as the mass ra- 
tio decreases. The decrease is expected since temperature 



is directly related to the mass of the particle (equations 
©). Thus, the lighter particles will have a lower temper- 
ature than the heavier particles. Component 2 exhibits 
a minimum as the height increases, just as in the case 
of component 1. However, the minimum becomes shal- 
lower as the temperature decreases. Figure^ also shows 
the temperature obtained from the packing fraction using 
Equation [5] Again, we find that there is good agreement 
between the temperature calculated by Equation [5] and 
that calculated from Equations for large heights. 



C. Effect of varying size ratio 

Finally, we studied the effects of particle size by sim- 
ulating systems with size ratios of d-^/dx = 1.0, 0.8, 
0.5, and 0.1. The mass ratio was held constant at 
mijm\ = 1.0, and the relative fraction was held con- 
stant at n^jn-i = 1.0, with n\ + ri2 = 1050. The changes 
observed in the energy distributions, the packing fraction 
and the partial temperatures are presented and discussed 
below. 



1. Energy Distribution 

Changes in particle size results in changes in the mean 
free path |27l and the pair correlation function at contact 
[Til ll5L llrJ]. These changes affect the particle velocities 
by changing the number of collisions that a particle ex- 
periences in a given amount of time. Thus, changes in 
the size ratio are expected to result in changes in the 
distribution of the energy, just as observed for the mass 
ratio. 

Figures ITUk and 1 10b show the average energy changes 
that particles of components 1 and 2 experience as a re- 
sult of collision. Collisions with the wall cause a loss 
of energy for both components while collision with base 
increase the energy of the particles for all size ratios (Fig- 
ur ell Obi. The collisions between particles, however show 
different trends than seen in the mass ratio study. As 
seen in Figure ITUb . the inter-component collisions do not 
result in an energy loss for component 1 for all size ra- 
tios. The smaller the size ratio, the greater the amount 
of energy injected per collision into component 1 from 
collisions with component 2. All particle-particle colli- 
sions decrease the energy of component 2. The energy 
loss per collision for for both inter-component and intra- 
component, particle-particle collisions increase as the size 
ratio decreases. It is interesting to note that the energy 
loss due to both inter-component and intra-component 
collisions are the same for component 2 at the smallest 
size ratio. The energy changes for all but the smallest size 
ratio agree well with those predicted by Equation[7| The 
smallest size ratio shows a large deviation between the en- 
ergy loss predicted by the equations and that determined 
from the simulation. One reason for the discrepancy is 
the assumption of Maxwell-Boltzmann velocities used in 
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FIG. 10: Average energy change per particle per collision 
for (a) particle-boundary collisions and (b) particle-particle 
collisions as the size ratio is changed. Component 1 is repre- 
sented by the closed symbols and component 2 is represented 
by the open symbols. The different kinds of collisions are: B 
intra-component, particle-particle ((Ai?i)i and (Ai?2)2), ^ 
inter-component, particle-particle ((AEi}^ and (A-E^i), ^ 
particle-wall, and A particle-base. The solid (dashed) lines 
are the energy loss for component 1 (component 2) as cal- 
culated by Equation for intra-component (same shade as 
the squares) and inter-component (same shade as the circles), 
particle-particle collisions. 

determining the equations. The differences in the dis- 
tribution of energy within each component and between 
the two components that we observe here will affect the 
packing fraction and temperature for these systems. 

2. Packing Fraction 

Figure ^2 shows the effect of the size ratio on the pack- 
ing fraction of the system. We can see that component 



1 reaches higher altitudes as the diameter of component 

2 decreases without any noticeable shift in the position 
of the maximum in packing fraction. This implies that 
the particles of component 1 expand through the system 
as the size ratio decreases. There are two possible causes 
of the increase in the tail of the packing fraction at large 
heights. First, the particles of component 1 are able to 
retain more energy because of a decrease in the collision 
rate between particles as the size ratio decreases. The 
decrease in the collision rates is the result of a decrease 
in the total excluded volume of the system as the size of 
component 2 is reduced (deduced from Figures ITTk and 
HTb). The change in the total excluded volume of the sys- 
tem reduces the amount of energy lost to particle-particle 
collisions. The other cause of the increased altitude is the 
energy gain that comes from collisions with component 2. 
Since energy can be gained at positions above the base, 
the particles will be able to travel to higher altitudes in 
the system. 

Figure HTb shows the effect of the size ratio on compo- 
nent 2. We can see an overall reduction in the packing 
fraction of component 2 as the particle size is reduced. 
This occurs because the number of particles of compo- 
nent 2 is held fixed as the size is decreased. We also 
see that the particles of component 2 are able to reach 
greater altitudes in the system as the size ratio decreases 
(see the insert in FigurelTTb). The increase in the altitude 
is the result of the decrease in the number of collisions 
discussed above. 



3. Temperature 

FigurefHlk shows the z-component of the granular tem- 
perature as a function of height for component 1 for the 
four size ratios. Figure 112b shows the same for compo- 
nent 2 of the mixture. A minimum is again observed 
in the temperatures of each system for each component. 
The other general features of the temperature profiles in 
the x- and y-directions, while not shown here, are the 
same as seen and discussed in connection with Figure 
El It can be easily seen in Figure El that the temper- 
ature of both the large and small particles increase as 
the size ratio decreases. This increase in temperature 
results from a decrease in the energy lost due to particle- 
particle collisions. The reduction in the collision rate as 
the size ratio decreases can be observed in Table H] It is 
interesting that the energy loss due to intra-component 
collisions increases for component 2 as the collision rate 
for intra-component collisions decreases (Figure ITUb ). In 
addition, there is a discrepancy between the theory and 
the simulation results observed in Figure ^| at small size 
ratios that should be noted. The deviations are the re- 
sult of the systems being dominated by collisions with 
the wall and not particle-particle collisions (see Table [IJ), 
as assumed in the theory (see the Appendix). 

We used Equation Inland the high altitude tails of the 
packing fractions shown in the inserts in Figures ITTk and 
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TABLE I: Collision rates (s 1 ) for the different kinds of collisions for components 1 and 2. 





Component 1 


Component 2 


Size Ratio 


intra-component 


inter-component 


Wall 


Base 


intra-component 


inter-component 


Wall 


Base 


1.0 


272.8 


546.2 


181.2 


78.8 


272.1 


546.2 


195.1 


78.8 


0.8 


252.4 


412.0 


278.5 


74.7 


164.1 


412.0 


381.2 


77.1 


0.5 


230.0 


248.1 


503.0 


70.3 


55.2 


248.1 


2587.4 


71.9 


0.1 


208.9 


72.5 


649.4 


68.1 


7.75 


72.5 


136270 


54.4 



0.04- 



0.04- 



0.02- 




100 150 
Height (mm) 



FIG. 11: The packing fraction of (a) 
component 2 for size ratios of B dz/di 

♦ da/di = 0.5, and A fc/di = 0.1. 
parameters are as given in the text. 



component 1 and (b) 
= 1.0, • da/di = 0.8, 
All other simulation 
The inset shows the 



natural log of the trailing edge of the packing fraction for 
each system. 



II lb to calculate a temperature. The temperatures calcu- 
lated in this manner are also plotted in Figure^] The 
temperature calculated from Equation El coincides with 
the temperature in the z-direction at large heights, just 
as in the case of the mass ratio. The agreement between 
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FIG. 12: The temperature in the z-direction of (a) component 
1 and (b) component 2 for size ratios of B di/d\ = 1.0, 
• da/dx = 0.8, ♦ d 2 /di = 0.5, and A fa/di = 0.1. The 

temperature calculated from Equation |5] is shown as a line 
corresponding to the data points of the same shade. 



the two temperatures is very good for both components 
for each size ratio, except the smallest. This may indicate 
that the packing fraction is not a very sensitive measure 
of the temperature when the difference in size between 
the components is large. The difference in the two tem- 
peratures may also be indicative of the changes in the 
dominant processes in the energy distribution. For exam- 
ple, the smallest size ratio is not dominated by particle- 
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particle collisions like the other systems. The edge effects 
associated with the walls of the cylinder become more im- 
portant, but they are not considered in the theory used 
to obtain Equation |5J 



IV. CONCLUSIONS 

We have shown that, like the single component system 
|18| , the simulation results reproduce the phenomena ob- 
served in experimental studies. Specifically, the experi- 
mentally observed flow pattern, radial dependence of 
the packing fraction and temperature of the two compo- 
nents are qualitatively reproduced by the simulation. 

There are a number of possible reasons for the quan- 
titative differences between the simulation and experi- 
ments. For example, for each component the simulated 
particles are identical spheres, while in the experimental 
system there is a small distribution of shape, size, and 
mass. The simulation also assumes that the particles are 
frictionless with constant restitution coefficients, which 
is not the case for ballotini glass spheres used in the ex- 
periments. The saw-tooth wave form of the vibrating 
base is an idealization of the sinusoidal vibration of the 
experiment. It is thought, however, that this does not 
have a large influence on the system behavior [3l]]. In 
any case, we feel confident that our model captures the 
key physical aspects of the experimental system and can 
therefore be used to study the influence of various system 
parameters. 

We also studied the effects of mass and size ratio in 
this paper. Generally, we observed that changes in either 
ratio do not result in any segregation of the particles. The 
lighter particles attain greater heights than the heavier 
particles in the mass ratio studies. For both components 
the rate of decay of the packing fraction at large heights 
is a good indication of the granular temperature in that 
region HHU- 

As the mass ratio decreases, the overall temperature 
of the system decreases. This is consistent with the 
lower amount of energy gained by the lighter compo- 
nent from the vibrating base and an overall lowering of 
the efficiency of energy transfer as the mass ratio de- 
creases. The packing fraction and temperature of the 
individual components, however, have a non-trivial de- 
pendence on the mass ratio. For the heavier compo- 
nent, both these quantities exhibit a minimum as the 
mass ratio decreases. We also observed that the en- 
ergy changes due to inter-component collisions exhibit 
extrema for each component as the mass ratio decreases. 
Specifically, the heavier component shows a minimum, 
whereas the lighter component shows a maximum in the 
energy change. The lighter component actually gains en- 
ergy from collisions with the other component around a 
mass ratio of m^/mi = 0.5. The energy changes due 
to intra-component collisions for both components, how- 
ever, are negative for all mass ratios. The energy loss due 
to intra-particle collisions increases for the heavy compo- 



nent as the mass ratio decreases, while it decreases for 
the lighter component. For comparison, we developed 
an approximate theory to calculate the energy changes 
for particle-particle collisions by assuming that the par- 
ticle velocities follow a Maxwcll-Boltzmann distribution 
with a component specific temperature. While this as- 
sumption is not strictly correct given the inhomogeneities 
in density and temperature, the theory is qualitatively 
accurate for the inter-component collisions and in near 
quantitative agreement with the simulation for the intra- 
component collisions. 

As the size ratio decreases at constant mass ratio, the 
overall temperature of the system increases and the to- 
tal packing fraction decreases. The particle-particle and 
particle-boundary collision rates decrease and increase, 
respectively, as the size ratio decreased. At the same 
time, the larger particles begin to gain energy from col- 
lisions with the smaller component. The approximate 
theory of the energy changes is again able to reproduce 
qualitatively the observed trends. However, the agree- 
ment between the theory and simulation results worsens 
as the size ratio decreases. 

Wildman and Parker observed a decrease in the 
temperature as they decreased the ratio of the number 
of large to small ballotini spheres. Our results show that 
this effect is dominated by the difference in mass of the 
two components, and not the difference in size. 
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APPENDIX: ENERGY DISSIPATION 
EQUATIONS 

The energy change of particle i of component a result- 
ing from collision with particle j of component (3 is 



2m a AE a 



fj 2 (1 + c) 2 (vy • ft) 2 



m a fi (1 + c) 



(v Q ,i • ft) 2 



foj • ft) 2 + ( Vij • ft) 2 (A.l) 



where vy = v a i — v^j and /z = m a mp/(m a +mp). The 
total energy loss due to the collision is simply AE to tai — 
AE a ,i + AE 0d , or 

A£ totQ; = -i/i(l-c 2 )(v ir fi) 2 (A.2) 

These equations are exact for each collision that occurs 
in the system. In order to determine the average values, 
however, some assumptions must be made. Specifically, 
we assume that the velocities of each component are de- 
scribed by a Maxwell-Boltzmann distribution that is ho- 



13 



mogeneous, isotropic and characterized by a component- 
specific temperature T a : 



f (v Q ,i) dv a ,i 



2itm a 
k B T a 



3/2 



exp 



2k B T a 



(A.3) 

In order to average over the velocities v Qj j, v^j, and vy 
appearing in Eg uat ions l A . 1 1 and l A . 21 we introduce center- 
of-mass, v c , and relative v r velocities: 



V Q ,i 



V /3J 



m a /T a + mp/Tp 

m a /T a ^ ^ 

m a/T a + 171(3 /Tp 

Vij = v Q: i - v/3j = v r 



(A.4) 

(A.5) 
(A.6) 



In this coordinate system the energy change for particle 
i becomes 

2m a AE a>i = ^ 2 (l + c) 2 (v r -n) 2 

—2T a fj,'[i (1 + c) (v r • n) 2 (A.7) 
-2m a [i (1 + c) (v r • n) (v c • n) 



where 



z _ {m a /T a ){mf3/Tf3) 



(A. 



m a /T a + m l3 /T l 3 
The total energy loss (shown in Eauation lA.2|) becomes 



A£ totai = -- A1 (l-c 2 )(v r -n) 2 (A.9) 



in the new coordinate system. 

It is now necessary to average Equations I A. 71 and IA.9I 

over the fraction of collisions between a particle of com- 
ponent a with a particle of component (3 with a relative 
velocity between v r and v r + dv r . Straightforward mod- 
ification of the standard kinetic theory result [221 gives: 



p( Vr )dv r = U-Q „3 exp ^_|g )(/r < A - IU1 



We then compute 



((v r -n) 2 }=/ dv r p(v r ) I dbh(b)(v r -hf (A.ll) 



where h(b)db = (2b/a 2 )db is the probability that the im- 
pact parameter lies between b and b + db. Substituting 
(v r • n) 2 = v 2 (l — (b/r) 2 ) and evaluating the integrals 
leads to 



it ~\ 2 \ 
v r • n = — 



(A.12) 



Using this result and the fact that ((v r • n) (v c • n)} =0 
gives Equations and |H1 
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